function CalcOpt()
    S = load('workParams.mat');
    
    t0 = 10000;
    ting = 365;        

    pStr2Double = defParamSet(S.pdm);

    resParams0 = CalcResParams(S.pf, pStr2Double, 0);    
    optsol0 = NNsysExact_ing0_p(t0, resParams0);
    statConc = deval(optsol0, t0);

    resParams15 = CalcResParams(S.pf, pStr2Double, 1.5e+6);    
    optsol15 = NNsysExact_ingNonzero_p(ting, resParams15, statConc);

    resParams75 = CalcResParams(S.pf, pStr2Double, 7.5e+6);    
    optsol75 = NNsysExact_ingNonzero_p(ting, resParams75, statConc);

    resParams150 = CalcResParams(S.pf, pStr2Double, 15e+6);    
    optsol150 = NNsysExact_ingNonzero_p(ting, resParams150, statConc);
    
    save('optsol.mat', 'optsol0', 'optsol15', 'optsol75', 'optsol150', 'statConc'); 
end